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ABSTRACT 

We describe a conservative, shock-capturing scheme for evolving the equations of 
general relativistic magnetohydrodynamics. The fluxes are calculated using the Harten, 
Lax, and van Leer scheme. A variant of constrained transport, proposed earlier by 
Toth, is used to maintain a divergence free magnetic field. Only the covariant form of 
the metric in a coordinate basis is required to specify the geometry. We describe code 
performance on a full suite of test problems in both special and general relativity. On 
smooth flows we show that it converges at second order. We conclude by showing some 
results from the evolution of a magnetized torus near a rotating black hole. 

Subject headings: accretion, accretion disks, black hole physics, Magnetohydrodynam- 
ics: MHD, Methods: Numerical 

1. Introduction 

Quasars, active galactic nuclei (AGN), X-ray binaries, gamma-ray bursts, and core-collapse su- 
pernovae are all likely powered by a central engine subject to strong gravity, strong electromagnetic 
fields, and rotation. For convenience, we will refer to this class of objects as relativistic magneto- 
rotators (RMRs). RMRs are among the most luminous objects in the universe and are therefore the 
center of considerable theoretical attention. Unfortunately the governing physical laws for RMRs, 
while well known, are nonlinear, time-dependent, and intrinsically multidimensional. This has 
stymied development of a first-principles theory for their evolution and observational appearance 
and strongly motivates a numerical approach. 
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To fully understand RMR structure, one must be able to follow the interaction of a non- 
Maxwellian plasma with a relativistic gravitational field, a strong electromagnetic field, and possibly 
a strong radiation field as well. This general problem remains beyond the reach of today's algorithms 
and computers. A useful first step, however, might be to study these objects in a nonradiating 
magnetohydro dynamic (MHD) model. In this case the plasma can be treated as a fluid, greatly 
reducing the number of degrees of freedom, and the radiation field can be ignored. The relevance 
of this approximation must be evaluated in astrophysical context and will not be considered here. 

We were motivated, therefore, to develop a method for integrating the equations of ideal, 
general relativistic MHD (GRMHD), and that method is described in this paper. This is in a sense 
well-trodden ground: many schemes have already been developed for relativistic fluid dynamics. 
What has not existed until recently is a scheme that: (1) includes magnetic fields; (2) has been 
fully verified and convergence tested; (3) is stable and capable of integrating a fiow over many 
dynamical times. A pioneering GRMHD code has been developed by Koide and collaborators (e.g. 
Koide et al. (2002), Koide, Shibata, & Kudoh (1999), Koide et al. (2000), Meier, Koide, & Uchida 
(2001)). Our code code differs in that we have subjected it to a fuller series of tests (described in 
this paper), we can perform longer integrations than the rather brief simulations described in the 
published work of Koide's group, and our code explicitly maintains the divergence-free constraint 
on the magnetic field. A GRMHD Godunov scheme based on a Roe-type approximate Riemann 
solver has been developed by Komissarov and described in a conference proceeding (Komissarov 
2001). 

A rather complete review of numerical approaches to relativistic fluid dynamics is given by 
Marti & Muller (1999) and Font (2000). The flrst numerical GRMHD scheme that we are aware 
of is by Wilson (1977), who integrated the GRMHD equations in axisymmetry near a Kerr black 
hole. While it was recognized throughout the 1970s that relativistic MHD would be relevant to 
problems related to black hole accretion (particularly following the seminal work of Blandford 
& Znajek (1977) and Phinney (1983); see also Punsly (2001)), no further work appeared until 
Yokosawa (1993). The discovery by Balbus &; Hawley (1991) that magnetic flelds play a crucial 
role in regulating accretion disk evolution (reviewed in Balbus k, Hawley (1998)) and the absence 
of purely hydrodynamic means for driving accretion in disks (Balbus & Hawley 1998) further 
motivated the development of relativistic MHD schemes. More recently there have been several 
efforts to develop GRMHD codes, including the already-mentioned work by Koide and collaborators 
and by Komissarov. A ZEUS-like scheme for GRMHD has also been developed and is described in 
a companion paper (de Villiers & Hawley 2002). 

Special relativistic MHD (SRMHD) is the foundation for any GRMHD scheme, although there 
are nontrivial problems in making the transition to full general relativity. SRMHD schemes have 
been developed by Van Putten (1993); Balsara (2001); Koldoba et al. (2002); Komissarov (1999) 
and Del Zanna et al. (2002). We were particularly influenced by the clear development of the 
fundamental equations in Komissarov (1999) for his Godunov SRMHD scheme based on a Roe- 
type approximate Riemann solver, and by the work of Del Zanna &l Bucciantini (2002) and Del 
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Zanna et al. (2002) who chose to use the simple approximate Riemann solver of Harten et al. (1983) 
in their special relativistic hydrodynamics and SRMHD schemes, respectively. 

Our numerical scheme is called HARM, for High Accuracy Relativistic Magnetohydrodynam- 
ics.^ In the next section we develop the basic equations in the form used for numerical integration 
in HARM (§2). In §3 we describe the basic algorithm. In §4 wc describe the performance of the 
code on a series of test problems. In §5 we describe a sample evolution of a magnetized torus near 
a rotating black hole. 



2. A GRMHD Primer 

The equations of general relativistic MHD are well known, but for clarity we will develop them 
here in the same form used in numerical integration. Unless otherwise noted c = 1 and we follow 
the notational conventions of Misner, Thorne, &; Wheeler (1973), hereafter MTW. The reader may 
also find it useful to consult Anile (1989). 

The first governing equation describes the conservation of particle number: 

[nu^^),^ = 0. (1) 

Here n is the particle number density and is the four- velocity. For numerical purposes wc rewrite 
this in a coordinate basis, replacing n with the "rest-mass density" p = mn, where m is the mean 
rest-mass per particle: 

-±=d^{V^pu>^) = 0. (2) 

Here g = Det{gi^t,). 

The next four equations express conservation of energy-momentum: 

T^;/. = 0, (3) 
where T'*!^ is the stress-energy tensor. In a coordinate basis, 

dt {V^ T\) = -di ( V=9 T\) + V=5 r\«, (4) 

where i denotes a spatial index and V^vk is the connection. 

The energy-momentum equations have been written with the free index down for a reason. 
Symmetries of the metric give rise to conserved currents. In the Kerr metric, for example, the 
axisymmetry and stationary nature of the metric give rise to conserved angular momentum and 
energy currents. In general, for metrics with an ignorable coordinate the source term on the 
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right hand side of the evolution equation for vanish. These source terms do not vanish when 
the equation is written with both indices up. 

The stress-energy tensor for a system containing only a perfect fluid and an electromagnetic 
field is the sum of a fluid part, 

TZd = iP + n+v)u>'u^ (5) 
(here u = internal energy and p = pressure) and an electromagnetic part, 

T^;^ = F'^^F^-\g'^'^F^f,F^^. (6) 

Here F^^'^ is the electromagnetic field tensor (MTW: "Faraday"), and for convenience we have 
absorbed a factor of -v/ivr into the definition of F. 

The electromagnetic portion of the stress-energy tensor simplifies if wc adopt the ideal MHD 
approximation, in which the electric field vanishes in the fluid rest frame due to the high conductivity 
of the plasma ( "E + v x B = 0" ) . Equivalently the Lorentz force on a charged particle vanishes in 
the fluid frame: 

u^F^'' = 0. (7) 
It is convenient to define the magnetic field four-vector 

where e is the Levi-Civita tensor. Recall that (following the notation of MTW) e'^'^^^ = — [ij.i/X6] , 
where [/xz^A(5] is the completely antisymmetric symbol and = 0, 1, or —1. These can be combined 
(with the aid of identity 3.50h of MTW): 

= e^""^\jx- (9) 

Substitution and some manipulation (using identities 3.50 of MTW and b^uf^ = 0; the latter follows 
from the definition of b^ and the antisymmetry of F) yields 

T^^ = b'^u^u'' + ^b'^g"" - b^b''. (10) 

Notice that the last two terms are nearly identical to the nonrelativistic MHD stress tensor, while 
the first term is higher order inv/c. To sum up, 

^MHD = ip + u + p + b')u^'u^ + (p+ ]f)g^^^ - Ifb'' (11) 

is the MHD stress-energy tensor. 

The electromagnetic field evolution is given by the source-free part of Maxwell's equations 

Fij.v,\ + -pA/i,i^ + Fv\,ii. = 0. (12) 
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The rest of Maxwell's equations determine the current 

= F^^., (13) 

and are not needed for the evolution, as in nonrelativistic MHD. 

Maxwell's equations can be written in conservative form by taking the dual of eq.(12): 

F*'^'^;^ = 0. (14) 

Here F*^ = ^e^i^^A-f"*'^ is the dual of the electromagnetic field tensor (MTW: "Maxwell"). In ideal 
MHD 

F*^^ = h'^u'' - b'^u'', (15) 

which can be proved by taking the dual of eq.(9). 

The components ofb^ are not independent, since b^U/j, = 0. Following, e.g., Komissarov (1999), 
it is useful to define the magnetic field three- vector = F***. In terms of B^, 

6* = B'u^^Qi^, (16) 

b' = {B' + b*u')/uK (17) 
The space components of the induction equation then reduce to 

dtiV^B^) = -d,{^ {Vu' - 6V)) (18) 

and the time component reduces to 

1 



=0^(7^ i?') = 0, (19) 



V-5 

which is the no-monopoles constraint. The appearance of B^ in these last two equations are what 
motivates the introduction of the field three-vector in the first place. 

To sum up, the fundamental equations as used in HARM are: the particle number conservation 
equation (2); the four energy- momentum equations (4), written in a coordinate basis and using 
the MHD stress-energy tensor of equation (11); and the induction equation (18), subject to the 
constraint (19). These hyperbolic ^ equations are written in conservation form, and so can be 
solved numerically by well-known techniques. 



®The GRMHD equations exhibit the same degeneracies as the nonrelativistic MHD equations. 
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3. 



Numerical Scheme 



There are many possible ways to numerically integrate the GRMHD equations. A first, zeroth- 
order choice is between conservative and nonconservative schemes. Nonconservative schemes such 
as ZEUS (Stone & Norman 1992) have enjoyed wide use in numerical astrophysics. They permit 
the integration of an internal energy equation rather than a total energy equation. This can be 
advantageous in regions of a flow where the internal energy is small compared to the total energy 
(highly supersonic flows), which is a common situation in astrophysics. A nonconservative scheme 
for GRMHD following a ZEUS-like approach has been developed and is described in a companion 
paper (de Villiers & Hawley 2002). 

We have decided to write a conservative scheme. One advantage of this choice is that in one 
dimension, total variation stable schemes are guaranteed to converge to a weak solution of the 
equations by the Lax-Wendroff theorem (Lax k. Wendroff 1960) and by a theorem due to LeVeque 
(1998). While no such guarantee is available for multidimensional flows, this is a reassuring starting 
point. Furthermore, one is guaranteed that a conservative scheme in any number of dimensions 
will satisfy the jump conditions at discontinuities. This is not true in artificial viscosity based 
nonconservative schemes, which are also known to have trouble in relativistic shocks (Norman & 
Winkler 1986). Conservative schemes for GRMHD have also been developed by Komissarov (2001) 
and by Koide, Shibata, & Kudoh (1999). 

A conservative scheme updates a set of "conserved" variables at each timestep. Our vector of 
conserved variables is 



These are updated using fluxes F. We must also choose a set of "primitive" variables, which are 
interpolated to model the flow within zones. We use variables with a simple physical interpretation: 



Here = /u^ is the 3-velocity.^ The functions U(P) and F(P) are analytic, but the inverse 
operations (so far as we can determine) are not. ^ There is also no simple expression for F(U). 

To evaluate U(P) and F(P) one must find n* and h^^ from and B^. To find u*, solve the 
quadratic equation g^i^u^vy = —1. Next use equations (16) and (17) to find b^] these require only 
multiplications and additions. The remainder of the calculation of U(P) and F(P) requires raising 
and lowering of indices followed by direct substitution in equation (11) to find the components of 
the MHD stress-energy tensor. 



®We initially used as primitive variables, but the inversion is not always single- valued, e.g. inside the 

ergosphere of a black hole. That is, there are physical flows with the same value of but different values of u*. 

^Del Zanna et al. (2002) have found that the inversion P(U) can be reduced analytically to the solution of a 
single, nonlinear equation. 



\5 = ^{pu\TlTlB'). 



(20) 



V = {p,u,v\B'). 



(21) 
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Since we update U rather than P, we must solve for P(U) at the end of each timestep. We 
use a multidimensional Newton-Raphson routine with the value of P from the last timestep as 
an initial guess. Since can be obtained analytically, only 5 equations need to be solved. The 
Newton-Raphson method requires an expensive evaluation of the Jacobian d\J/dP. In practice we 
evaluate the Jacobian analytically. It is possible to evaluate the Jacobian by numerical derivatives, 
but this is both expensive and a source of numerical noise. 

The evaluation of P(U) is at the heart of our numerical scheme; the procedure must be robust. 
We have found that it is crucial that the errors (differences between the current and target values 
of U) used to evaluate convergence in the Newton-Raphson scheme be properly normalized. We 
normalized the errors with ^—gpu*. 

To evaluate F we use a MUSCL type scheme with "HLL" fluxes (Harten et al. 1983). The 

fluxes are defined at zone faces. A slope-limited linear extrapolation from the zone center gives 
P/j and Pi, the primitive variables at the right and left side of each zone interface. Wc have 
implemented the monotonized central ("Woodward", or "MC") limiter, the van Leer limiter, and 
the minmod limiter; unless otherwise stated, the tests described here use the MC limiter, which is 
the least diffusive of the three. 

From P/j, P^, calculate the maximum left and rightgoing wave speeds c±^r, c±^l, and the fluxes 
Fr = F(Pr) and Fl = F(Pi). Defining Cmax = MAX(0, c+,r, c+,l) and c^i„ = -MIN(0, c_,r, c_,l), 
the HLL flux is then 

■pi (^min^R "I" ^max^ L CmaxCmini^ R (22) 

Cmax "l~ (Wain 

If Cmax = CmiH) the HLL flux becomes the so-called local Lax-Priedrichs flux. 



3.1. Constrained Transport 

The pure HLL scheme will not preserve any numerical representation of V • B = 0. An incom- 
plete list of options for handling this constraint numerically includes: (1) ignore the production of 
monopoles by truncation error and hope for the best (in our experience this causes the scheme to 
fail in any complex flow); (2) introduce a divergence-cleaning step (this entails solving an elliptic 
equation at each timestep); (3) use an Evans & Hawley (1988) type constrained transport scheme 
(this requires a staggered mesh, so that the magnetic held components are zone face centered); (4) 
introduce a diffusion term that causes numerically generated monopoles to diffuse away (Harder 
(1987); this typically leaves a monopole field with rms value somewhat larger than the truncation 
error) . 

We have chosen a fifth option, a version of constrained transport that can be used with a 
zone-centered scheme. This idea was introduced by one of us in Toth (2000), where it is called 
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the flux- interpolated constrained transport (or "flux-CT") scheme ^. It preserves a numerical 
representation of V • B = by smoothing the fluxes with a special operator. The disadvantage of 
this method is that it is more diffusive than the "bare", unconstrained scheme. The advantage is 
that it is extremely simple. 

To clarify how zone-centered constrained transport works we now give a specific example for 
a special relativistic problem in Cartesian coordinates t, x, y. To fix notation, write the induction 
equation as 

dtB' = -djF^, (23) 
where we have used = 1 and the fluxes Ff are 



^ X 

Fl 

px 

y 






yuyx 



l)X^y 

l,yyx 



(24) 



-Fl 



Notice that the fluxes are centered at different locations on the grid: F^ fluxes live on the x face of 
each zone at grid location i — 1/2, j (we use i,j to denote the center of each zone), while F^ fluxes 
live on the y face at i,j — 1/2. The smoothing operator replaces the numerical (HLL-derived) F- 
with , deflned by 



4"(i-V2,j) 

F^hj-m 
FUij-m 








(25) 



l[2Flii,j-l/2) 

+Fl{i + 1, j - 1/2) + Fl{i - 1, j - 1/2) 
-F-{i-l/2,j)-F-{i + l/2,j) 
-F-{i - 1/2, j - 1) - F-{i + 1/2, j - 1) 

F-{i- 1/2, j) = l[2F-ii-l/2,j) 

+F^{i - 1/2, j + 1) + F^{i - 1/2, j - 1) 
-F|(i,j-l/2)-F|(i,j + l/2) 
-F|(i -l,j- 1/2) - F|(i -l,j + 1/2) 

It is a straightforward but tedious exercise to verify that this preserves the following corner-centered 
numerical representation of V • B: 



V B 



(5-(i,i) + B-(i,i - 1) - B-ii - l,j) - B-{i -1,3- 1)) /(2Ax) 



+ {By{i,j) + By{i - l,j) - By{i,j - l) - By{i -l,j- l)) /(2Ay), 

where A,T and Ay arc the grid spacing. 



(26) 



We have also experimented with Toth's "flux-CD" , which preserves a different representation of V • B = 0. This 
appears to be shghtly less robust. It also has a larger effective stencil. 
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3.2. Wave Speeds 



The HLL approximate Riemann solver does not require eigenvectors of the characteristic ma- 
trix (as would a Roe-type scheme), but it does require the maximum and minimum wave speed 
(eigenvalues). These wave speeds are also required to fix the timestep via the Courant conditions. 

The relevant speed is the phase speed "w/A;" of the wave, and it turns out that only speeds for 
waves with wavevectors aligned along coordinate axes are required. Suppose, for example, that 
one needs to know how rapidly signals propagate in the fluid along the xi direction. First, find 
a wavevector k/^ = (— a;, fei, 0, 0), that satisfies the dispersion relation for the mode in question: 
D{k^) = 0. Then the wave speed is simply oo/ki. 

The dispersion relation D^k^) = for MHD waves has a simple form in a comoving frame. 
In terms of the relativistic sound speed = {d{p + u)/dp)j^ = jp/w, (the last holds only if 
p = (7— 1)-^) and the relativistic Alfven velocity va = B/\/J, where 6 = b^ + w and w = p + u+p, 
the dispersion relation is 



Here c is the (temporarily reintroduced) speed of light. The first term is the zero frequency entropy 
mode, the second is the Alfven mode, and the third contains the fast and slow modes. The eighth 
mode is eliminated by the no-monopoles constraint. 

The relativistic sound speed asymptotes to c-\/j — 1 = c/\/3 for 7 = 4/3, and the Alfven 
speed asymptotes to c. In the limit that B'^/p » 1 and p/p 1, the GRMHD equations 

are a superset of the time-dependent, force-free electrodynamics equations recently discussed by 
Komissarov (2002b); these contain fast modes and Alfven modes that move with the speed of 
light. They are indistinguishable from vacuum electromagnetic modes only when their wavevector 
is oriented along the magnetic field. 

To find the maximum wave speeds we need to evaluate the comoving-frame dispersion relation 
for the fast wave branch from coordinate frame quantities. This is straightforward because the 
dispersion relation depends on scalars, which can be evaluated in any frame: uj = kfj_u^; k^ = K^K^, 
where = {g^v+u^Ui,)k^ is the part of the wavevector normal to the fluid 4-velocity; = h^J)^ jS; 
(k • va) = kp^^ l^fE. The relevant portion of the dispersion relation (for fast and slow modes) is 
thus a fourth order polynomial in the components of k^. This can be solved either analytically or 
by standard numerical methods. The two fast mode speeds are then used as Cjjiax cind Cniin 

the 

HLL fluxes. 

We have found it convenient to replace the full dispersion relation by an approximation: 



This overestimates the maximum wavespeed by a factor < 2 in the comoving frame. The maximum 
error occurs for k || va,va/cs = 1, and va <^ c, and it is usually much less, particularly if the fluid 



u (w^ — (k ■ va)^) X 

{u^ + (A;2(v2 + c2(l - vi/c^)) + cj{k ■ va)^^') + k^cjO^ ■ va)^) = 0, 



(27) 



u;' = ivl + cl{l-vl/c'))k\ 



(28) 
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is moving super-Alfvenically with respect to the grid. This approximation is convenient because it 
is quadratic in k^, and so can be solved more easily. 

3.3. Implementation Notes 

For completeness we now give some details of the implementation of the algorithm. 

Time Stepping. Our scheme is made second order in time by taking a half-step from to 
evaluating F(P(t"+V2))^ and using that to update U(t") to U(t"+^). 

Modification of Energy Equation. A direct implementation of the energy equation can be 

inaccurate because the magnetic and internal energy density can be orders of magnitude smaller 
than the rest mass density. To avoid this we subtract the particle number conservation equation 
from the energy equation, i.e., we evolve 

dtiV^in + pn*)) = -di{^g{Ti + pu')) + V^gnvl. (29) 

In the nonrelativistic limit, this procedure subtracts the rest mass energy density from the total 
energy density. 

Specification of Geometric Quantities. In two dimensions we need to evaluate g^u,9^'^, 
and ^/—g at four points in every grid zone (the center, two faces, and one corner) and F^^ at the 
zone center. It would be difficult to accurately encode analytic expressions for all these quantities. 
HARM is coded so that an analytic expression need only be provided for g^,^; all other geomet- 
ric quantities are calculated numerically. The connection, for example, is obtained to sufficient 
accuracy by numerical differentiation of the metric. This minimizes the risk of coding errors in 
specifying the geometry. It also minimizes coordinate dependent code, making it relatively easy to 
change coordinate systems. Minimal coordinate dependence, besides following the spirit of general 
relativity, enables one to perform a sort of fixed mesh refinement by adapting the coordinates to 
the problem at hand. For example, near a Kerr black hole we use log(r) as the radial coordi- 
nate instead of the usual Boyer-Lindquist r, and this concentrates numerical resolution toward the 
horizon, where it is needed. 

Density and Internal Energy Floors. Negative densities and internal energies are for- 
bidden by the GRMHD equations, but numerically nothing prevents their appearance. In fact, 
negative internal energies are common in numerical integrations with large density or pressure con- 
trast. Following common practice, we prevent this by introducing "floor" values for the density 
and internal energy. These floors are enforced after the half-step and the full step. They preserve 
velocity but do not conserve rest mass or energy-momentum. 

Outflow Boundary Conditions. In the rotating black hole calculations described below we 
use outflow boundary conditions at the inner and outer radial boundaries. The usual implementa- 
tion of outflow boundary conditions is to simply copy the primitive variables from the boundary 
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zones into the ghost zones. This can result in unphysical values of the primitive variables in the 
ghost zones- for example, velocities that lie outside the light cone- because of variations in the 
metric between the boundary and ghost zones. 

We have experimented with a variety of schemes for projecting variables into the ghost zones 
in the context of black hole accretion flow calculations (described in §5). We find that some are 
more robust than others. The most robust extrapolates the density, internal energy, and radial 
magnetic field according to 

P(ghost) = P(boundary)-v/^(boundary)/\/^(ghost), (30) 

the 9 and (j) components of the velocity and magnetic field according to 

P(ghost) = P(boundary)(l. - Ar/r), (31) 
and the radial velocity according to 

P(ghost) = P(boundary)(l. + Ar/r). (32) 

The extrapolation of 9 and (j) components of magnetic field and velocity results in weak damping 
of these components near the boundary. Slightly different choices of the extrapolation coefficients 
(i.e. (1. — 2Ar/r)) are much less robust. 

Performance. We have implemented both serial and parallel versions of the code. In serial 
mode the code integrates the black hole accretion problem (described in §5) at « 54, 000 zone cycles 
per second on a 2.4 GHz Intel Pentium 4, when compiled using the Intel C compiler. The parallel 
code was implemented using MPI. 

4. Code Verification 

Here we present a test suite for verifying a GRMHD code. The tests are nonrelativistic, special 
and general relativistic, and one and two dimensional. The list of problems for which there are 
known, exact solution is short, since exact solutions of multidimensional GRMHD problems are 
algebraically complicated. This list of test problems was developed in collaboration with J. Hawley 
and J.-P. de Villiers. Unless otherwise stated we set 7 = 4/3 and c = 1. 

4.1. Linear Modes 

This first test considers the evolution of a small amplitude wave in two dimensions. The 
unperturbed state is p = 1, p = 1, = 0, 5^ = 5^ = 0, 5^ = . The basic state is parametrized 
by a = (Bq)^/(/9c^); our fiducial test runs have a = 1. Onto this basic state we introduce a 
perturbation of the form exp(zk • x — ia;(k)t), where {kx,ky) = (27r, 27r), and the amplitude is fixed 
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by dB^ = 10 ^-Bq. The computational domain is x,y e [0, 1), [0, 1), and the boundary conditions 
are periodic. The wave is either slow, Alfvenic, or fast. 

This test exercises almost all terms in the governing equations. The numerical resolution is 

(Nx,Ny) = N{5,4:) zones, and the integration runs for a single wave period 27r/a;, so that a perfect 
scheme would return the simulation to its original state. We measure the Ci norm of the difference 
between the final state and the initial state for each primitive variable. For example, we measure 



for the density. 

All primitive variables exhibit similar convergence properties (as they must, since with the 
exception of the magnetic field, they are tightly coupled together). In Figures 1, 2, and 3 we 
present the Ci norm of the error for runs using the monotonized central limiter and a Courant 
number of 0.8, in addition to the results for the minmod limiter. These runs have a = 1. Figure 1 
shows the results for the slow wave. Figure 2 for the Alfven wave, and Figure 3 for the fast wave. 
Evidently the convergence rate asymptotes to second order, although more slowly for the minmod 
limiter. 

The code performs similarly well over a range of a, provided only that 5B^/2 <^ p, which is 
necessary for the wave to be in the linear regime. We have been unable to find a value of Bq where 
the code fails completely for a linear amplitude disturbance, although for very large values of Bq 
the evolution becomes inaccurate because of numerical noise in the evaluation of P(U). 



Komissarov (1999) has proposed a suite of one-dimensional nonlinear tests for special relativis- 
tic MHD. Komissarov presents a total of 9 tests (see his Table 1). The nonlinear Alfven wave (test 
5), and the compound wave (test 6) cannot be reconstructed without a separate derivation of the 
exact analytic solution, and we will not provide that here. For the remaining tests Komissarov's 
Table 1 contains several misprints that are corrected in Komissarov (2002a). Our code is able 
to integrate each of Komissarov's remaining 7 tests, although in some cases we must reduce the 
Courant number (usually 0.8) or resort to the slightly more robust van Leer slope limiter. Tests 
that required special treatment are: fast shock (Courant number = 0.5); shock tube 1 (Courant 
number = 0.3, van Leer limiter); shock tube 2 (Courant number = 0.5); collision (Courant number 
= 0.3; van Leer limiter). 

Figures 4 and 5 show the run of p and u^, respectively, for all 7 tests. These may be compared 
with Komissarov's figures. Notice that, unlike Komissarov, we have in all cases set = 400 and 
X G (—2, 2). We have not obtained the exact solutions used by Komissarov, but the solutions can 
still be checked quantitatively. For example, the slow shock wave speed is 0.5; since the calculation 
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Nonlinear Waves 
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ends at t = 2 the slow shock front should be, and is, located at x ^ 1.0. The fast shock speed is 
0.2, so at t = 2.5 the fast shock wave front should be, and is, located dX x ^ 0.5. 

There are artifacts evident in the figures. In particular there is ringing near the base of the 
switch-on and switch-off rarefaction waves. This is common and is seen in Komissarov's results as 
well. In addition the narrow, Lorentz-contracted shell of material behind the shock in shock tube 
1 is poorly resolved; the correct shell density is « 0.88 but a resolution of > 800 is required to 
find this result to an accuracy of a few percent. There is also a transient associated with the fast 
shock that propagates off the grid and so is not visible in Figures 4 and 5. 

A complete set of nonlinear wave tests for one dimensional nonrelativistic MHD was developed 
by Ryu & Jones (1995) (hereafter RJ). We can run these under HARM by rescaling the speed of 
light to c = 10^ in code units, where all velocities in the tests are 0{1). This should lead to results 
that agree with RJ to 0{v/c) 1%. The results can be checked quantitatively by comparison to 
the tables provided by RJ. 

Figure 6 shows our results for RJ test 5A, which is a version of the familiar Brio & Wu (1988) 
magnetized shock tube test. Like RJ we use 512 zones between ,t = and x = 1, and wc measure 
the results at t = 0.15. To compare to RJ quantitatively, consider Ux in the region behind the 
fast rarefaction wave, near x = 0.7. RJ report Ux = —0.277 here, while we measure Ux = —0.273, 
which differs by 1%, as expected. Similar agreement is found for the other variables. The most 
unsatisfactory feature of the solution is the visible post-shock oscillations. The amplitude of these 
features varies, depending on the Courant number (here 0.9) and the choice of slope limiter (here 
MC). 

Figure 7 shows our results for RJ test 2A. In the region near x = 0.6, RJ report that By = 
1.4126, and we find essentially exact agreement {By = 1.41262) after averaging over a small region 
near x = 0.6. This test has two pairs of closely spaced slow shocks and rotational discontinuities 
that are difficult to resolve, and our scheme barely obtains the correct peak values of Uy and By, 
even though there are about 20 zones inside the "horns" visible in the By panel of the figure. 

4.3. Transport 

This special relativistic test evolves a disk of enhanced density moving at an angle to the 
grid until it returns to its original position. The computation is carried out in a domain x,y E 
[—0.5, 0.5), [—0.5, 0.5) and the boundary conditions are periodic. The initial state has = = 0.7, 
or = ^ 4.95, corresponding to u* ~ 7.07. The initial density p = 1 except in a disk at 
r < Vg = 0.45, where p = 3/2 -|- cos(27rr/rs). The initial pressure p = 1, and the initial magnetic 
field is zero. The test is run until t = 10/7, when the system should return exactly to its initial 
state. 

Numerically, we use the monotonized central limiter and set the Courant number to 0.8. The 
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resolution is fixed so that N^, = 5Ny/A. Figure 8 shows the Li norm of the error in p as a function 
of X resolution. The convergence rate asymptotes to second order. 



The Orszag-Tang vortex (OTV) is a classic nonlinear MHD problem (Orszag &: Tang 1979). 
Here we compare our code, with the speed of light set to 100 so that that it is effectively non- 
relativistic, to the output of VAC (Toth 1996), an independent nonrelativistic code developed by 
one of us. The version of VAC used here is TVD-MUSCL using the monotonized central limiter. 
It is dimensionally unsplit and uses a scheme similar to HARM to control V • B. The problem is 
integrated in the periodic domain x G (— tt, tt], y E (— vt, tt] from t = to t = tt. Our version of the 
OTV has 7 = 4/3, but is otherwise identical to the standard problem. 

Results are shown in Figure 9, which shows p along a cut through the model at y = 7r/2 and 
t = TT. The resolution is 640^. The solid line shows the results from HARM; the dashed line shows 
the results from VAC. The lower solid line shows the difference between the two multiplied by 4. 
Evidently our code behaves similarly to VAC on this problem. 

We can quantify this by asking how the difference between the HARM and VAC solutions 
changes as a function of resolution. Figure 10 shows the variation in the Ci norm of the difference 
between the two solutions. Thus the line marked p shows 



evaluated at i = tt. The codes converge to one another approximately linearly, as expected for a 
flow containing discontinuities. If this study were extended to higher resolution convergence would 
eventually cease because the HARM solution would differ from the VAC solution due to finite 
relativistic corrections. 



Spherically symmetric accretion (Bondi flow) in the Schwarzschild geometry has an analytic 
solution (see, e.g., Shapiro &; Teukolsky (1983)) that can be compared with the output of our code. 
This appears to be a one-dimensional test, but for HARM it is actually two dimensional. Although 
the pressure is independent of the Boyer-Lindquist coordinate 9, the 9 acceleration does not vanish 
identically. This is because pressure enters the momentum equations through a flux {—de{psm9) in 
the Newtonian limit) and a source term (p cos 9 in the Newtonian limit) . Analytically these terms 
cancel; numerically they produce an acceleration that is of order the truncation error. 

Our test problem follows that set out in Hawley, Smarr, &; Wilson (1984): we fix the sonic 
point Ts = SGM/c^, M = Airr^pu'^ = —1, and 7 = 4/3. The problem is integrated in the domain 



4.4. Orszag-Tang Vortex 
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4.5. Bondi Flow in Schwarzschild Geometry 
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r e (1.9,20)GM/c^ for At = lOOGM/c^. We use coordinates based on the Kerr-Schild system, 
whose hne element is 

ds^ = -(l-2r/p2)(it2^(4r/p2)a!rdt + (l + 2r/p2)dr2 + pW+ 

sin^ (p2 + a2 ( 1 + 2r /p2 ) sin^ 0) #2 (35) 
-(4ar sin^ e/p'^)dtd(p - 2a{l + 2r / p^) sin^ edrdcj), 

where we have set GM = c = 1. In (35) only = + a^cos^(0); elsewhere p is density. In this 

test, a = 0. Wc modify these coordinates by replacing r by xi = log(r). The new coordinates are 
implemented by changing the metric rather than changing the spacing of grid zones. We measure 
the £1 norm of the difference between the initial conditions (exact analytic solution) and the final 
state. The difference is taken over the inner 3/4 of the grid in each direction, thus excluding 
boundary zones where errors may scale differently. This test exercises many terms in the code 
because in Kerr-Schild coordinates only three of the ten independent components of the metric are 
zero. 

The C\ norm of the error in internal energy for the Bondi test is shown in Figure 11. Similar 
results obtain for the other independent variables. The solution converges at second order. 

4.6. Magnetized Bondi Flow 

The next test considers a Bondi flow containing a spherically symmetric, radial magnetic field. 
The solution to this problem is identical to the Bondi fiow described above because the flow is 
along the magnetic field, so all magnetic forces cancel exactly. This is a difficult test, however, 
because numerically the magnetic terms cancel only to truncation error. This causes problems at 
high magnetic field strength. 

We use the same Bondi solution as in the last subsection and parameterize the magnetic field 
strength by at the inner boundary. Our fiducial run has {b^/p){r = l.QGM/c^) = 10.56. The 
Ci norm of the error in the internal energy is shown in Figure 12. Similar results obtain for the 
other independent variables. The solution converges at second order. 

We have considered models with a range of (6^/p)(ri„). Lowering (6^/p)(ri„) produces results 
similar to those in our fiducial test run. Raising {b^/p){rin) first causes the code to produce 
inaccurate results (at ~ 10^, where the radial velocity profile is smoothly distorted from the true 
solution) and then to fail (at ~ 10^). This is an example of a general problem with conservative 
schemes when the basic energy density scales (rest mass, magnetic, and internal) differ by many 
orders of magnitude. 
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4.7. Magnetized Equatorial Inflow in Kerr Geometry 

This test considers the steady-state, magnetized inflow solutions found by Takahashi et al. 
(1990), as specialized to the case of inflow inside the marginally stable orbit by Gammie (1999). 
This solution exercises many of the important terms in the governing equations, in particular the 
interaction of the magnetized fluid with the Kerr geometry. 

We use Boyer-Lindquist coordinates to specify this problem, but the solution is integrated in the 
modified Kerr-Schild coordinates described above. The flow is assumed to lie in the neighborhood 
of the black hole's equatorial plane and is thus one dimensional, much like the Weber-Davis model 
for the solar wind. As above, we set GM = c = 1. 

The particular inflow solution we consider is for a black hole with spin parameter a/M = 0.5. 
The model has an accretion rate Fm = — 1 = 27r pr'^u'^ (adopting the notation of Gammie 1999). 
The magnetization parameter Fg^ = r^B^ = 0.5. The flow is constrained to match to a circular 
orbit at the marginally stable orbit. This is enough to uniquely specify the flow. It follows that 
(see Gammie 1999) Ftg = ClFg^, where O is the orbital frequency at the marginally stable orbit. 
For a/M = 0.5, CI 0.10859. The solution that is regular at the fast point has angular momentum 
flux Fl = 2'Kr'^{vJ'u^ - Ifh^) « -2.8153 and energy flux Fe = 2'Kr'^{vrut - Wht) « -0.90838. 
The fast point is located at r ~ 3.6167, and the radial component of the four-velocity there is 
= —0.040547. Figure 13 shows the radial run of the solution. 

We initialize the flow with a numerical solution that is subject to roundoff error. The near- 
equatorial nature of the solution is mimicked by using a single zone in the direction centered on 
^ = 7r/2. The computational domain runs from 1.02x the horizon radius r/j to 0.98x the radius of 
the marginally stable orbit Vmso- For a/M = 0.5, = 1.866, and rmso = 4.233. The analytic flow 
model is cold (zero temperature) but we set the initial internal energy in the code equal to a small 
value instead. The model is run for At = 15. 

Figure 14 shows the jCi norm of the error in p, u'^,u'^, and B'^ and a function of the total number 
of radial gridpoints N. The straight line shows the slope expected for second order convergence. 
The small deviation from second order convergence at high N in several of the variables is due to 
numerical errors in the initial solution, which relies on numerical derivatives (Gammie 1999). 

4.8. Equilibrium Torus 

Our next test concerns an equilibrium torus. This class of equilibria, found originally by 
Fishbone & Moncrief (1976) and Abramowicz, Jaroszinski, &: Sikora (1978), consist of a "donut" of 
plasma surrounding a black hole. The donut is supported by both centrifugal forces and pressure 
and is embedded in a vacuum. Here we consider a particular instance of the Fishbone & Moncrief 
solution. 
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A practical problem with this test is that HARM abhors a vacuum. We have therefore intro- 
duced floors on the density and internal energy that limit how small these quantities can be. The 
floors are dependent on radius, with Pmin = 10^^(r/rj„)~^/^ and Umin = 10~^(r/rj„)~^/^. This 
means that the torus is surrounded by an insubstantial, but dynamic, accreting atmosphere that 
interacts with the torus surface. To minimize the influence of the atmosphere on our convergence 
test, we take the Ci norm of the change in variables only over that region where p > Q.Q2pTnax- 

The problem is integrated in modified Kerr-Schild coordinates. The Kerr-Schild radius r has 
been replaced by the logarithmic radial coordinate xi = ln(r), and the Kerr-Schild latitude 9 
has been replaced by X2 such that 6 = 7rx2 + (1/2)(1 — /i) sin(27rx2). Clearly < X2 < 1 maps 
to < < TT. This coordinate transformation has a single adjustable parameter h] for = 1 
we recover the original coordinate system (the 9 coordinate is simply rescaled by tt). As /i — > 
numerical resolution is concentrated near the midplane. 

Wc have integrated a Fishbonc-Moncricf disk around a black hole with a/M = 0.95, to max- 
imize general relativistic effects. Wc set u^Ufj) = const. = 3.85 (this is the defining feature of the 
Fishbone-Moncrief equilibria) and Tin = 3.7. The grid extends radially from Rin = 0.98r/i = to 
Rout = 20. The coordinate parameter h described in the last paragraph is set to 0.2. The numerical 
resolution is NxN, where N = 8, 16, 32, ... , 512, and the solution is integrated for At = 10. Figure 
15 shows the Ci norm of the error for each variable as a function of N. Second order convergence 
is obtained. 

The sum of the evidence presented in this section strongly suggests that we are solving the 
equations of GRMHD without significant, compromising errors. 



5. Magnetized Torus Near Rotating Black Hole 

Finally we offer an example of how HARM can be applied to a real astrophysical problem: the 
evolution of a magnetized torus near a rotating black hole. Again wc set GM = c = 1. 

The initial conditions contain a Fishbone-Moncrief torus with a/M = 0.5, r{pmax) = 12, and 
Tin = 6. Superposed on this equilibrium is a purely poloidal magnetic field with vector potential 
OC MAX(p/ Pmax ~ 

0.2,0), where 

Pmax is the peak density in the torus. The field is normalized so 
that the minimum value of Pgas/Pmag = 10^- The orbital period at the pressure maximum (r = 12), 
is 264 as measured by an observer at infinity. 

The integration extends for At = 2000, or about 7.6 orbital periods at the pressure maximum. 
Figure 16 shows the initial and final density states projected on the R = rsm{9), Z = rcos{9) 
plane. Color represents log{p). The coordinate parameter h, which concentrates zones toward the 
midplane, is set to 0.2. The torus atmosphere is set to the floor values (see above), and the MC 
limiter is used. The numerical resolution is 300^. 

The flux of mass, energy, and angular momentum through the inner boundary are described 
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in Figure 17. Initially the fluxes are small because the initial conditions are near an (unstable) 
equilibrium. The magnetorotational instability (Balbus &: Hawley 1991) e- folds for just over an 
orbital period, after which the magnetic field has reached sufficient strength to distort the original 
torus and drop material into the black hole. Later, the torus is turbulent and accretion occurs at 
a more or less steady rate. 

6. Conclusion 

Like all hydrodynamics codes, HARM has failure modes. We will discuss one that is likely 
to be relevant to future astrophysical simulations. When -B^/p S> 1 and 5^ S> u, the magnetic 
energy is the dominant term in the total energy equation. Because the fields arc evolved separately, 
truncation error in the field evolution can lead to large fractional errors in the velocity and internal 
energy. An example of this was discussed in §4.6, where the magnetized Bondi flow test fails for 
large values of B^/p. 

Another example can be found in the strong cylindrical explosion problem of Komissarov 
(1999), where an overpressured region embedded in a uniform magnetic field produces a relativistic 
blast wave. HARM fails on the strong-field version of this problem unless we turn the Courant 
number down to 0.1, use the minmod limiter, and sharply increase the accuracy parameter used in 
the P(U) inverter. This is a particularly difficult problem, with B'^/p as large as 10'^. The problems 
caused by magnetically dominated regions appears to be generic to conservative relativistic MHD 
schemes, where small errors in magnetic energy density lead to fractionally large errors in other 
components of the total energy. At present this is unavoidable, and has motivated the development 
of schemes for the evolution of the electromagnetic field in the force-free limit (Komissarov 2002b) . 

Finally, to sum up: we have described and tested a code that evolves the equations of general 
relativistic magnetohydrodynamics. This code, together with the code described in a companion 
paper by de Villiers & Hawley (2002) , are the first that stably evolve a relativistic plasma in a Kerr 
spacetime for many light crossing times. The advent of practical, stable GRMHD codes opens the 
door for the study of many problems in the theory of RMRs. For example, it may be possible to 
directly evaluate the importance of magnetic energy extraction from rotating black holes and the 
importance of black hole spin in determining jet parameters. It may also be possible to couple 
these schemes to numerical relativity codes and use them to study dynamical spacetimes with 
electromagnetic sources. 

We are grateful to our collaborators John Hawley and Jean-Pierre de Villiers for extensive 
discussions. We are also grateful to Ramesli Narayan, who supplied the initial inspiration for this 
project some years ago. Stu Shapiro's advice has greatly improved this paper. Serguei Komissarov 
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Fig. 1. — The Ci norm of the error in u for a slow wave as a function of A^^ for both the monotonized 
central (MC) and minmod limiter. The straight lines show the slope expected for second order 
convergence. 
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Fig. 2. — The L\ norm of the error in the single nonzero component of the velocity for an Alfven 
wave as a function of for both the monotonizcd central (MC) and minmod limiter. The straight 
lines show the slope expected for second order convergence. 
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Fig. 3. — The Ci norm of the error in u for a fast wave as a function of for both the monotonized 
central (MC) and minmod hmiter. The straight Unas show the slope expected for second order 
convergence. 
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The run of density in the Komissarov nonlinear wave tests. 
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Fig. 5.— 



The run of in the Komissarov nonhnear wave tests. 
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Fig. 6. — Snapshot of the final state in HARM's integration of Ryu &; Jones test 5A (a version 
of the Brio & Wu shock tube) but with c = 100. The figure shows primitive variable values at 
t = 0.15. Quantitative agreement is found to within 1%, as expected. 
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Fig. 7. — Snapshot of the final state in HARM's integration of Ryu & Jones test 2A, with c = 100. 
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Convergence results for the transport test. 
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X (at y = 7r/2, i = tx) 

Fig. 9. — A cut through the density in the nonrelativistic Orszag-Tang vortex solution from HARM 
(solid line, with c = 100), from VAC (dashed line), and 4x the difference (lower solid line) at a 
resolution of 640^. 
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Fig. 10. — Comparison of results from HARM and the nonrelativistic MHD code VAC for the 

Orszag-Tang vortex. The plot shows the C\ norm of the difference between the two results as a 
function of resolution for the primitive variables p (squares) and u (triangles). The straight line 
shows the slope expected for first order convergence. The errors are large because they are an 
integral over an area of (27r)^. 
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Fig. 11. — Convergence results for the unmagnetized Bondi accretion test onto a Schwarzschild 
black hole. The straight line shows the slope expected for second order convergence. 
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Fig. 12. — Convergence results for the magnetized Bondi accretion test onto a Schwarzschild black 
hole. The straight line shows the slope expected for second order convergence. 
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Fig. 13. — The equatorial inflow solution in the Kerr metric for a/M = 0.5 and magnetization 
parameter Fg^ = 0.5. The panels show density, radial component of the four-velocity in Boyer- 
Lindquist coordinates (with the square showing the location of the fast point), the ^ component of 
the four- velocity, and the toroidal magnetic field = F'l'*. 
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Fig. 14. — Convergence results for the magnetized inflow solution in a Kerr metric with a/M = 0.5. 
Parameters for the initial, quasi-analytic solution are given in the text. The straight line shows 
the slope expected for second order convergence. The Ci error norm for each of the nontrivial 
variables are shown. The small deviation from second order convergence at high resolution is due 
to numerical errors in the quasi-analytic solution used to initialize the solution. 
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Fig. 15. — Convergence results for the Fishbone and Moncrief equihbrium disk around an a/M = 
0.95 black hole. 
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Fig. 16. — Density field, for a magnetized torus around a Kerr black hole with a/M = 0.5 at t = 
(left) and at t = 2000M (right). The color is mapped from the logarithm of the density; black is 
low and dark red is high. The resolution is 300^. 
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Fig. 17. — Evolution of the mass accretion rate (top), the specific energy of the accreted matter 
(middle) , and the specific angular momentum of the accreted matter (bottom) for a black hole with 
a/M = 0.5. 



